{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 量子变分求解器在量子化学计算中的应用\n",
    "\n",
    "\n",
    "`Linux` `CPU` `全流程` `初级` `中级` `高级`\n",
    "\n",
    "[![](https://gitee.com/mindspore/docs/raw/master/tutorials/training/source_zh_cn/_static/logo_source.png)](https://gitee.com/mindspore/mindquantum/blob/master/tutorials/vqe_for_quantum_chemistry.ipynb)\n",
    "\n",
    "## 概述\n",
    "\n",
    "量子化学，指的是运用量子力学的基本理论及方法，求解含时或定态薛定谔方程的数值解。在高性能计算机上进行量子化学模拟已成为研究材料的物理、化学性质的重要手段。然而，精确求解薛定谔方程具有指数级的复杂度，可模拟的化学体系规模严重受制于此。近年量子计算的发展为解决这个问题提供了一条可行的路，有望在量子计算机上实现多项式复杂度下对薛定谔方程的高精度求解。\n",
    "\n",
    "[Peruzzo等人](https://doi.org/10.1038/ncomms5213)在2014年首次将量子变分求解器(Variational quantum eigensolver, VQE)结合[幺正耦合簇理论](https://doi.org/10.1016/S0009-2614(89)87372-5)用于量子化学的模拟中，实现了He-H<sup>+</sup>基态能量的求解。量子变分求解器是一个量子-经典混合算法，在基于量子算法的化学模拟中应用广泛，本教程将介绍使用量子变分求解器求解分子体系基态能量的方法。\n",
    "\n",
    "本教程的主要内容包括如下几个部分：\n",
    "\n",
    "1. 量子化学原理简介\n",
    "2. 量子变分求解器的应用\n",
    "3. 使用mindquantum实现高效自动求导的VQE模拟\n",
    "\n",
    "## 环境准备\n",
    "\n",
    "本教程需要安装以下环境：\n",
    "- NumPy\n",
    "- SciPy\n",
    "- [mindquantum](https://gitee.com/mindspore/mindquantum)\n",
    "- [mindspore](https://gitee.com/mindspore/mindspore)\n",
    "- PySCF\n",
    "- openfermion\n",
    "- openfermionpyscf\n",
    "\n",
    "导入本教程所依赖模块"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86fbb794",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from openfermion.chem import MolecularData\n",
    "from openfermionpyscf import run_pyscf\n",
    "import mindquantum as mq\n",
    "from mindquantum import Circuit, X, Hamiltonian\n",
    "from mindquantum.circuit import generate_uccsd\n",
    "from mindquantum.nn import generate_pqc_operator\n",
    "import mindspore as ms\n",
    "import mindspore.context as context\n",
    "from mindspore.common.parameter import Parameter\n",
    "from mindspore.common.initializer import initializer\n",
    "\n",
    "context.set_context(mode=context.GRAPH_MODE, device_target=\"CPU\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2f89dfd",
   "metadata": {},
   "source": [
    "## 量子化学计算方法\n",
    "\n",
    "量子化学的核心问题在于求解薛定谔方程(Schrödinger Equation)。一般来说，求解含时薛定谔方程(Time-dependent Schrödinger Equation)较为复杂，故引入玻恩-奥本海默近似(Born-Oppenheimer approximation, BO approximation)。BO近似认为，原子核质量远大于电子、运动速度远低于电子，故可以将两者进行分离变量，单独讨论原子核或电子的运动，于是可得到如下不含时的电子运动方程，也称为定态薛定谔方程：\n",
    "\n",
    "$$\n",
    "\\hat{H} |\\Psi\\rangle = E |\\Psi\\rangle\n",
    "$$\n",
    "\n",
    "其中$\\hat{H}$包含以下三项：\n",
    "\n",
    "$$\n",
    "\\hat{H} = \\hat{K} _{e} + \\hat{V} _{ee} + \\hat{V} _{Ne}\n",
    "$$\n",
    "\n",
    "分别为电子动能、电子-电子势能和电子-核势能。\n",
    "\n",
    "有多种数值算法可以求解定态薛定谔方程。本教程将介绍其中的一类：波函数方法。波函数方法直接求解给定分子哈密顿量的本征波函数和本征能量，目前有大量的开源软件包可实现，如[PySCF](http://pyscf.org/)等。此处从一个简单的例子：氢化锂分子开始，使用openfermion结合openfermionpyscf插件进行。首先定义分子结构："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "1ef0fb6d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Geometry: \n",
      " [['Li', [0.0, 0.0, 0.0]], ['H', [0.0, 0.0, 1.5]]]\n"
     ]
    }
   ],
   "source": [
    "dist = 1.5\n",
    "geometry = [\n",
    "    [\"Li\", [0.0, 0.0, 0.0 * dist]],\n",
    "    [\"H\",  [0.0, 0.0, 1.0 * dist]],\n",
    "]\n",
    "basis = \"sto3g\"\n",
    "spin = 0\n",
    "print(\"Geometry: \\n\", geometry)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a6071dd",
   "metadata": {},
   "source": [
    "上面的代码定义了一个Li-H键长为1.5Å分子。使用STO-3G基组进行计算。接下来使用openfermionpyscf，调用PySCF进行HF、CCSD和FCI计算。这三种方法属于波函数方法，开始计算之前，先对这些方法作一个简单的介绍：\n",
    "\n",
    "### 波函数方法\n",
    "\n",
    "求解定态薛定谔方程的方法之一是[Hartree-Fock(HF)](https://doi.org/10.1098/rspa.1935.0085)方法，该方法在二十世纪三十年代左右由Hartree等人提出，是量子化学计算中的基本方法。HF方法引入了单行列式近似，即$N$-电子体系的波函数由一个行列式形式的波函数表示：\n",
    "\n",
    "$$\n",
    "| \\Psi \\rangle = | \\psi_{1} \\psi_{2} \\psi_{3} \\dots \\psi_{N} \\rangle\n",
    "$$\n",
    "\n",
    "其中$| \\psi_{1} \\psi_{2} \\psi_{3} \\dots \\rangle$代表由一组自旋轨道波函数$\\{ \\pi_{i} \\}$构成的N阶行列式。\n",
    "自旋轨道波函数$\\psi_{i}$可进一步用一组形式已知的基函数展开：\n",
    "\n",
    "$$\\psi_{i} = \\phi_{i} \\eta_{i}$$\n",
    "$$\\phi_{i} = \\sum_{\\mu}{C_{\\mu i} \\chi_{\\mu}}$$\n",
    "\n",
    "其中$\\{\\chi_{\\mu}\\}$被称为基函数，可以是高斯函数等。\n",
    "该近似考虑了电子间的交换作用，但是忽略了电子间的关联作用，故在无法正确计算如解离能等性质。\n",
    "\n",
    "HF方法的改进可以从波函数展开定理出发。波函数展开定理可以表述为，若$\\{ \\psi_{i} \\}$是一组完备的自旋轨道波函数，则$N$-电子体系波函数可以由$\\{ \\psi_{i} \\}$构成的行列式波函数精确展开：\n",
    "\n",
    "$$\n",
    "| \\Psi \\rangle = \\sum^{\\infty} _ {i_{1} < i_{2} < \\dots < i_{N}} {C_{i_{1} i_{2} \\dots i_{N}} | \\psi_{i_{1}} \\psi_{i_{2}} \\dots \\psi_{i_{N}} \\rangle}\n",
    "$$\n",
    "\n",
    "由此可得到Configuration Interaction(CI)方法：\n",
    "\n",
    "$$\n",
    "| \\Psi_{CI} \\rangle = C_{0} | \\Psi_{HF} \\rangle + \\sum^{a\\rightarrow\\infty} _{i\\in occ\\\\\\\\a\\not\\in occ}{C^{a} _{i} | \\Psi^{a} _{i} \\rangle } + \\sum^{ab\\rightarrow\\infty} _{ij\\in occ\\\\\\\\ab\\not\\in occ}{C^{ab} _{ij} | \\Psi^{ab} _{ij} \\rangle } \n",
    "$$\n",
    "\n",
    "上式中的$| \\Psi^{a}_{i} \\rangle + \\dots$代表电子由轨道$i$激发到轨道$a$的单激发波函数，以此类推。只考虑单激发和双激发的CI被称为CISD，即Configuration Interaction with singles and doubles。将基态HF波函数一直到N激发波函数全部考虑在内的Configuration Interaction被称为Full Configuration Interaction(FCI)，FCI波函数是定态薛定谔方程在给定基函数下的精确解。\n",
    "\n",
    "### 二次量子化\n",
    "\n",
    "在二次量子化表述下，体系的哈密顿量具有如下形式：\n",
    "\n",
    "$$\n",
    "\\hat{H} = \\sum_{p, q}{h^{p} _ {q} E^{p} _ {q}} + \\sum_{p, q, r, s}{\\frac{1}{2} g^{pq} _ {rs} E^{pq} _ {rs} }\n",
    "$$\n",
    "\n",
    "其中$E^{p} _ {q}$和$E^{pq} _ {rs}$分别为：\n",
    "\n",
    "$$\n",
    "E^{p} _ {q} = a^{\\dagger} _ {p} a_{q}\n",
    "$$\n",
    "$$\n",
    "E^{pq} _ {rs} = a^{\\dagger} _ {p} a^{\\dagger} _ {q} a_{r} a_{s}\n",
    "$$\n",
    "\n",
    "$a^{\\dagger} _ {p}$和$a _ {q}$分别为产生算符(Creation Operator)和湮灭算符(Annihilation Operator)。\n",
    "\n",
    "使用二次量子化的表述方法，可以非常方便地表示激发态波函数：\n",
    "\n",
    "$$\n",
    "| \\Psi^{abc\\dots} _ {ijk\\dots} \\rangle = a^{\\dagger} _ {a} a^{\\dagger} _ {b} a^{\\dagger} _ {c} \\dots a_{i} a_{j} a_{k} \\dots | \\Psi \\rangle\n",
    "$$\n",
    "\n",
    "CI方法的一个改进是耦合簇理论(Coupled-Cluster theory, CC)。CC引入指数化算符：\n",
    "\n",
    "$$\n",
    "| \\Psi_{CC} \\rangle = \\exp{(\\hat{T})} | \\Psi_{HF} \\rangle\n",
    "$$\n",
    "\n",
    "其中耦合簇算符$\\hat{T}$为对激发算符的求和：\n",
    "\n",
    "$$\n",
    "\\hat{T} = \\sum_{p\\not\\in occ\\\\\\\\q\\in occ}{\\theta^{p} _ {q} E^{p} _ {q}} + \\sum_{pq\\not\\in occ\\\\\\\\rs\\in occ}{\\theta^{pq} _ {rs} E^{pq} _ {rs}} + \\dots\n",
    "$$\n",
    "\n",
    "其中$\\theta$和CI方法中的$C$类似，是待求解的参数。由指数的泰勒展开易知，即使耦合簇算符$\\hat{T}$中只包含低阶激发项，$\\exp{(\\hat{T})}$也可以隐含部分高阶激发，这也使得CC方法向FCI波函数收敛的速度要远快于CI，同样截断到K激发，如K=2，CCSD的精度会超过CISD。\n",
    "\n",
    "\n",
    "<!--\n",
    "一般而言，若一个方法可以达到化学精度，即由此方法计算的能量和FCI能量之间的差值小于1 kcal/mol，则认为这个方法具有良好的精度，截断到三激发的CCSD(T)在大部分情况下都能符合这个标准\n",
    "-->\n",
    "\n",
    "电子关联作用的效果是使得总能量降低，故HF得到的基态能量会略高于CCSD和FCI。另外，从上述理论不难发现，FCI的计算量远大于CCSD和HF。我们使用openfermion封装的MolecularData和openfermionpyscf封装的run_pyscf函数来进行演示："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "302e3a7d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Hartree-Fock energy:  -7.8633576215351164 Ha\n",
      "CCSD energy:  -7.8823529091527016 Ha\n",
      "FCI energy:  -7.8823622867987249 Ha\n"
     ]
    }
   ],
   "source": [
    "molecule_of = MolecularData(\n",
    "    geometry,\n",
    "    basis,\n",
    "    multiplicity=2 * spin + 1\n",
    ")\n",
    "molecule_of = run_pyscf(\n",
    "    molecule_of,\n",
    "    run_scf=1,\n",
    "    run_ccsd=1,\n",
    "    run_fci=1\n",
    ")\n",
    "\n",
    "print(\"Hartree-Fock energy: %20.16f Ha\" % (molecule_of.hf_energy))\n",
    "print(\"CCSD energy: %20.16f Ha\" % (molecule_of.ccsd_energy))\n",
    "print(\"FCI energy: %20.16f Ha\" % (molecule_of.fci_energy))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8ef9ee4",
   "metadata": {},
   "source": [
    "在上面的例子中，我们运行了Hartree-Fock(HF)、CCSD、FCI进行总能量的计算。若对运行时间进行统计，会发现$T_{HF}<T_{CCSD}\\ll T_{FCI}$，换成计算量更大的体系如乙烯分子等会更明显一些。此外，对于计算得到的总能量，有$E_{HF}>E_{CCSD}>E_{FCI}$。计算完成后，我们将结果保存到`molecule_file`文件（即`molecule_of.filename`）中："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "633173b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "molecule_of.save()\n",
    "molecule_file = molecule_of.filename\n",
    "print(molecule_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd9a431f",
   "metadata": {},
   "source": [
    "量子化学计算的一大阻碍是计算量。随着体系大小（电子数、原子数）的增加，求解FCI波函数和基态能量的时间消耗大约以$2^{N}$增长，即使是较小的分子如乙烯分子等，进行FCI计算也并不容易。量子计算机的出现为此提供了一条可能的解决途径，已有的研究表明，量子计算机可以多项式的时间复杂度模拟哈密顿量的含时演化，在量子处理器上进行化学模拟相较于经典计算机有指数级的加速。本教程将介绍其中一类量子算法：量子变分求解器。\n",
    "\n",
    "## 量子变分求解器\n",
    "\n",
    "量子变分求解器(Variational Quantum Eigensolver, VQE)是一类量子-经典混合(Hybrid quantum-classical)算法，应用变分原理实现对基态波函数的求解。其中，变分参数的优化步在经典计算机上进行。\n",
    "\n",
    "### 变分原理\n",
    "\n",
    "变分原理可使用如下形式表述：\n",
    "\n",
    "$$\n",
    "E_{0} \\le \\frac{\\langle \\Psi_{t} | \\hat{H} | \\Psi_{t} \\rangle}{\\langle \\Psi_{t} | \\Psi_{t} \\rangle}\n",
    "$$\n",
    "\n",
    "上式中的$| \\Psi_{t} \\rangle$代表试探波函数。变分原理表明，在满足一定的条件下，任意试探波函数得到的基态能量总是大于等于真实的基态能量。变分原理为求解分子基态薛定谔方程提供了一种方法：使用一个参数化的函数$f(\\theta)$作为精确基态波函数的近似，通过优化参数$\\theta$来逼近精确的基态能量。\n",
    "\n",
    "### 初态制备\n",
    "\n",
    "在二次量子化表述下，$N$-电子HF波函数也具有非常简洁的形式：\n",
    "\n",
    "$$\n",
    "| \\Psi_{HF} \\rangle = \\prod^{i=0} _{N-1}{a^{\\dagger} _{i}| 0 \\rangle}\n",
    "$$\n",
    "\n",
    "上式搭建了一个由量子化学波函数到量子计算的桥梁：用$|0\\rangle$代表非占据轨道，用$|1\\rangle$代表电子占据的轨道，由此可以将$N$-电子HF波函数映射为由一串$M+N$个量子比特$| 00\\dots 11\\dots \\rangle$，$M$代表非占据轨道的数量。\n",
    "\n",
    "以下代码构造了对应于LiH分子的HF初态波函数。在Jordan-Wigner变换下，相当于将$N$个$\\text{X}$门作用于$|000\\dots\\rangle$上。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "087f484c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X(0)\n",
      "X(1)\n",
      "X(2)\n",
      "X(3)\n"
     ]
    }
   ],
   "source": [
    "hartreefock_wfn_circuit = Circuit([X.on(i) for i in range(molecule_of.n_electrons)])\n",
    "print(hartreefock_wfn_circuit)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef879afa",
   "metadata": {},
   "source": [
    "基于此，我们可以构造如下形式的试探波函数：\n",
    "\n",
    "$$\n",
    "| \\Psi_{t} \\rangle = U(\\theta) | \\Psi_{HF} \\rangle\n",
    "$$\n",
    "\n",
    "其中$U(\\theta)$代表一个可通过量子线路模拟的幺正变换，$| \\Psi_{HF} \\rangle$作为初态，可通过多个单比特$\\text{X}$门来方便地制备。$U(\\theta) | \\Psi_{HF} \\rangle$的具体形式也被称为波函数拟设。\n",
    "### 波函数拟设\n",
    "\n",
    "前文提到的耦合簇理论是一个非常高效的波函数拟设。在量子计算机上使用，需要作一些修改：\n",
    "\n",
    "$$\n",
    "| \\Psi_{UCC} \\rangle = \\exp{(\\hat{T} - \\hat{T}^{\\dagger})} | \\Psi_{HF} \\rangle\n",
    "$$\n",
    "\n",
    "UCC即幺正耦合簇(Unitary Coupled-Cluster theory)，$\\hat{T}^{\\dagger}$代表$\\hat{T}$的厄米共轭。如此，$\\exp{(\\hat{T} - \\hat{T}^{\\dagger})}$即为幺正算符。[Peruzzo等人](https://doi.org/10.1038/ncomms5213)在2014年首次使用VQE结合UCCSD(Unitary coupled-cluster with singles and doubles)拟设进行了量子计算机上的化学模拟实验。值得注意的是幺正耦合簇默认了耦合簇算符中的参数$\\{\\theta\\}$是实数。在分子体系中该假设不会有问题；在周期性体系中，[刘杰等人](https://doi.org/10.1021/acs.jctc.0c00881)的研究表明幺正耦合簇会因为忽略复数部分而造成误差。本教程暂时不讨论幺正耦合簇在周期性体系中的应用。\n",
    "\n",
    "使用mindquantum的circuit模块中的相应函数可读取先前保存在`molecule_file`的计算结果，方便地构造UCCSD波函数拟设，以及其对应的量子线路："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "0d14403e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ccsd:-7.882352909152702.\n",
      "fci:-7.882362286798725.\n"
     ]
    }
   ],
   "source": [
    "ansatz_circuit, \\\n",
    "init_amplitudes, \\\n",
    "ansatz_parameter_names, \\\n",
    "hamiltonian_QubitOp, \\\n",
    "n_qubits, n_electrons = generate_uccsd(molecule_file, th=-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1cd51664",
   "metadata": {},
   "source": [
    "完整的量子线路包含HF初态+UCCSD拟设，如下代码所示："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "43d02711",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===============================Circuit Summary===============================\n",
      "|Total number of gates  : 12612.                                            |\n",
      "|Parameter gates        : 640.                                              |\n",
      "|with 44 parameters are : p41, p28, p38, p10, p13, p31, p16, p3, p4, p29... |\n",
      "|Number qubit of circuit: 12                                                |\n",
      "=============================================================================\n",
      "Number of parameters: 44\n"
     ]
    }
   ],
   "source": [
    "total_circuit = hartreefock_wfn_circuit + ansatz_circuit\n",
    "total_circuit.summary()\n",
    "print(\"Number of parameters: %d\" % (len(ansatz_parameter_names)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21e8c28e",
   "metadata": {},
   "source": [
    "对于LiH分子而言，其UCCSD波函数拟设中包含44个变分参数。该线路总共的量子比特门数量为12612，总共需要12个量子比特进行模拟。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abaf223b",
   "metadata": {},
   "source": [
    "### VQE的一般流程\n",
    "\n",
    "使用VQE进行分子基态求解的一般流程如下：\n",
    "\n",
    "1. 制备HF初态：$| 00\\dots11\\dots \\rangle$；\n",
    "2. 定义波函数拟设，如UCCSD等；\n",
    "3. 将波函数拟设转化为参数化的量子线路；\n",
    "4. 初始化变分参数，如全设为0等；\n",
    "5. 在量子计算机上多次测量得到分子哈密顿量在该套变分参数下的能量$E(\\theta)$以及能量关于参数的导数$\\{ {\\partial E} / {\\partial \\theta_{i}} \\}$\n",
    "6. 在经典计算机上使用优化算法，如梯度下降、BFGS等更新变分参数；\n",
    "7. 将新的变分参数传入量子线路中进行更新；\n",
    "8. 重复步骤(5)到(7)，直到满足收敛标准；\n",
    "9. 结束\n",
    "\n",
    "在第5步中，求取能量关于参数的导数$\\{ {\\partial E} / {\\partial \\theta_{i}} \\}$在量子计算机上可通过parameter-shift rule来进行，在模拟器中也可通过模拟parameter-shift rule或者有限差分法来计算，是个较为耗时的过程。mindquantum基于mindspore框架，提供了类似于机器学习的自动求导功能，可以对参数化的量子线路进行反向传播求导，如此在模拟中可以达到高效计算导数的效果。以下使用mindquantum构造带自动求导功能的参数化UCCSD量子线路："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "6666fdb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "molecule_pqc = generate_pqc_operator(\n",
    "    [], ansatz_parameter_names, total_circuit, \n",
    "    Hamiltonian(hamiltonian_QubitOp))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39351f3e",
   "metadata": {},
   "source": [
    "通过将参数的具体数值传入`molecule_pqc`，即可得到对应于此变分参数的能量$E(\\theta)=\\langle \\Psi_{UCC}(\\theta) | \\hat{H} | \\Psi_{UCC}(\\theta) \\rangle$以及关于每个变分参数的导数。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87bd0e09",
   "metadata": {},
   "source": [
    "接下来需要进行VQE优化的(5)~(7)步，即对参数化量子线路进行优化。我们可以借助mindspore框架，使用参数化量子线路算子`molecule_pqc`构造一个机器学习模型，然后通过类似于训练神经网络的方法来优化变分参数。在[附录](#appendix)中会演示使用SciPy进行优化。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "cb2865a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "class PQCNet(ms.nn.Cell):\n",
    "    def __init__(self, pqc):\n",
    "        super(PQCNet, self).__init__()\n",
    "        self.pqc = pqc\n",
    "        self.weight =  Parameter(initializer(\"Zeros\",\n",
    "            len(self.pqc.ansatz_params_names)),\n",
    "            name=\"weight\")\n",
    "        self.encoder_data_dummy = ms.Tensor([[0]], \n",
    "            self.weight.dtype)\n",
    "\n",
    "    def construct(self):\n",
    "        energy, _, grads = self.pqc(self.encoder_data_dummy, self.weight)\n",
    "        return energy\n",
    "\n",
    "molecule_pqcnet = PQCNet(molecule_pqc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "此处使用`\"Zeros\"`关键字，将所有的变分参数初始化为0。使用CCSD（耦合簇理论）或者MP2（二阶多体微扰论）的计算结果也可以作为幺正耦合簇变分参数的初始值，此处暂不对此进行讨论。此时有$E(\\vec{0})=\\langle \\Psi_{UCC}(\\vec{0}) | \\hat{H} | \\Psi_{UCC}(\\vec{0}) \\rangle = E_{HF}$："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "0752a45c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial energy:  -7.8633575439453125\n"
     ]
    }
   ],
   "source": [
    "initial_energy = molecule_pqcnet()\n",
    "print(\"Initial energy: %20.16f\" % (initial_energy.asnumpy()))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6815cbef",
   "metadata": {},
   "source": [
    "最后使用mindspore的Adam优化器进行优化，学习率设置为$1\\times 10^{-2}$，优化终止标准设置为$\\left.|\\epsilon|\\right. = \\left.|E^{k+1} - E^{k}|\\right. \\le 1\\times 10^{-8}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "7411644c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Step   0 energy  -7.8633575439453125\n",
      "Step   1 energy  -7.8578791618347168\n",
      "Step   2 energy  -7.8738608360290527\n",
      "Step   3 energy  -7.8696637153625488\n",
      "Step   4 energy  -7.8716387748718262\n",
      "Step   5 energy  -7.8722548484802246\n",
      "Step   6 energy  -7.8763413429260254\n",
      "Step   7 energy  -7.8799986839294434\n",
      "Step   8 energy  -7.8786134719848633\n",
      "Step   9 energy  -7.8745274543762207\n",
      "Step  10 energy  -7.8739027976989746\n",
      "Step  11 energy  -7.8769011497497559\n",
      "Step  12 energy  -7.8790454864501953\n",
      "Step  13 energy  -7.8783979415893555\n",
      "Step  14 energy  -7.8772087097167969\n",
      "Step  15 energy  -7.8777532577514648\n",
      "Step  16 energy  -7.8795628547668457\n",
      "Step  17 energy  -7.8810067176818848\n",
      "Step  18 energy  -7.8813152313232422\n",
      "Step  19 energy  -7.8808994293212891\n",
      "Step  20 energy  -7.8806061744689941\n",
      "Step  21 energy  -7.8808226585388184\n",
      "Step  22 energy  -7.8812174797058105\n",
      "Step  23 energy  -7.8812923431396484\n",
      "Step  24 energy  -7.8810596466064453\n",
      "Step  25 energy  -7.8808951377868652\n",
      "Step  26 energy  -7.8809628486633301\n",
      "Step  27 energy  -7.8811511993408203\n",
      "Step  28 energy  -7.8813695907592773\n",
      "Step  29 energy  -7.8816132545471191\n",
      "Step  30 energy  -7.8818173408508301\n",
      "Step  31 energy  -7.8819088935852051\n",
      "Step  32 energy  -7.8819212913513184\n",
      "Step  33 energy  -7.8819446563720703\n",
      "Step  34 energy  -7.8820166587829590\n",
      "Step  35 energy  -7.8820791244506836\n",
      "Step  36 energy  -7.8820466995239258\n",
      "Step  37 energy  -7.8819408416748047\n",
      "Step  38 energy  -7.8819050788879395\n",
      "Step  39 energy  -7.8820013999938965\n",
      "Step  40 energy  -7.8821077346801758\n",
      "Step  41 energy  -7.8821139335632324\n",
      "Step  42 energy  -7.8820934295654297\n",
      "Step  43 energy  -7.8821563720703125\n",
      "Step  44 energy  -7.8822617530822754\n",
      "Step  45 energy  -7.8822956085205078\n",
      "Step  46 energy  -7.8822484016418457\n",
      "Step  47 energy  -7.8822150230407715\n",
      "Step  48 energy  -7.8822431564331055\n",
      "Step  49 energy  -7.8822717666625977\n",
      "Step  50 energy  -7.8822517395019531\n",
      "Step  51 energy  -7.8822278976440430\n",
      "Step  52 energy  -7.8822493553161621\n",
      "Step  53 energy  -7.8822855949401855\n",
      "Step  54 energy  -7.8822970390319824\n",
      "Step  55 energy  -7.8822999000549316\n",
      "Step  56 energy  -7.8823142051696777\n",
      "Step  57 energy  -7.8823270797729492\n",
      "Step  58 energy  -7.8823256492614746\n",
      "Step  59 energy  -7.8823184967041016\n",
      "Step  60 energy  -7.8823161125183105\n",
      "Step  61 energy  -7.8823189735412598\n",
      "Step  62 energy  -7.8823204040527344\n",
      "Step  63 energy  -7.8823213577270508\n",
      "Step  64 energy  -7.8823251724243164\n",
      "Step  65 energy  -7.8823308944702148\n",
      "Step  66 energy  -7.8823361396789551\n",
      "Step  67 energy  -7.8823409080505371\n",
      "Step  68 energy  -7.8823447227478027\n",
      "Step  69 energy  -7.8823432922363281\n",
      "Step  70 energy  -7.8823409080505371\n",
      "Step  71 energy  -7.8823409080505371\n",
      "Optimized energy:  -7.8823409080505371\n",
      "Optimized amplitudes: \n",
      " [ 3.7435064e-04  2.0003314e-03  3.4515299e-02  1.5759237e-02\n",
      "  2.2138040e-04  8.8742195e-04  5.9856582e-05  1.4393384e-02\n",
      "  6.4204258e-05  8.8728970e-04 -1.3264270e-04  1.4399976e-02\n",
      " -4.6512191e-04  5.3341745e-04  2.7783324e-03  5.2846380e-02\n",
      "  1.2344023e-04 -2.0551911e-04  5.8422484e-05  6.4405147e-05\n",
      " -5.9763537e-05 -2.5284235e-04 -2.4150473e-04 -5.0080267e-05\n",
      " -6.0626466e-05 -1.9624823e-04  1.3753232e-04  9.2733640e-04\n",
      " -5.4294106e-02  3.0179413e-03  6.6223278e-05 -7.7102435e-05\n",
      "  1.8220932e-04 -1.5686324e-04  6.5045329e-05  1.8605946e-06\n",
      "  1.6173313e-04  1.6479526e-04  3.0232230e-03 -8.6284184e-05\n",
      "  1.3787865e-04 -1.2463883e-04 -1.9412280e-04  4.1091337e-04]\n"
     ]
    }
   ],
   "source": [
    "optimizer = ms.nn.Adam(molecule_pqcnet.trainable_params(), learning_rate=1e-2)\n",
    "train_pqcnet = ms.nn.TrainOneStepCell(molecule_pqcnet, optimizer)\n",
    "\n",
    "eps = 1.e-8\n",
    "energy_diff = eps * 1000\n",
    "energy_last = initial_energy.asnumpy() + energy_diff\n",
    "iter_idx = 0\n",
    "while (abs(energy_diff) > eps):\n",
    "    energy_i = train_pqcnet().asnumpy()\n",
    "    print(\"Step %3d energy %20.16f\" % (iter_idx, float(energy_i)))\n",
    "    energy_diff = energy_last - energy_i\n",
    "    energy_last = energy_i\n",
    "    iter_idx += 1\n",
    "\n",
    "print(\"Optimized energy: %20.16f\" % (energy_i))\n",
    "print(\"Optimized amplitudes: \\n\", molecule_pqcnet.weight.asnumpy())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 附录\n",
    "\n",
    "<a id=\"appendix\"></a>\n",
    "\n",
    "目前流行的教程中，对于二次量子化的分子哈密顿量以及幺正耦合簇算符的具体细节介绍较为缺乏，给入门带来了一定的困难。本教程接下来的附录将对哈密顿量中的单、双电子积分$h^{p} _ {q}$和$g^{pq} _ {rs}$，以及耦合簇算符的具体计算和构造进行解释。以下将跳过openfermionpyscf，不使用其封装好的函数，而是直接从PySCF的计算结果出发，从头构造VQE计算所需要的所有物理量。\n",
    "\n",
    "接着前文的内容，导入一些额外的依赖："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "248d41b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from functools import reduce\n",
    "\n",
    "from scipy.sparse.linalg import eigsh\n",
    "from scipy.optimize import minimize\n",
    "import pyscf\n",
    "from pyscf import scf, cc, fci, ao2mo\n",
    "from mindquantum.ops import FermionOperator, QubitOperator\n",
    "from mindquantum.hiqfermion.transforms import Transform\n",
    "from mindquantum.utils import hermitian_conjugated, count_qubits, normal_ordered\n",
    "from mindquantum.circuit import decompose_single_term_time_evolution"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2a51dde",
   "metadata": {},
   "source": [
    "手动使用PySCF进行HF、CCSD和FCI计算："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "a8c72aab",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running Restricted Hartree-Fock...\n",
      "converged SCF energy = -7.86335762153512\n",
      "Running Restricted CCSD...\n",
      "E(CCSD) = -7.882352909152702  E_corr = -0.01899528761758478\n",
      "Running Full Configuration Interaction...\n",
      "Hartree-Fock energy:  -7.8633576215351164 Ha\n",
      "CCSD energy:  -7.8823529091527016 Ha\n",
      "FCI energy:  -7.8823622867987213 Ha\n"
     ]
    }
   ],
   "source": [
    "molecule = pyscf.gto.M(\n",
    "    atom=geometry,\n",
    "    basis=basis,\n",
    "    spin=spin\n",
    ")\n",
    "energy_nuc = molecule.energy_nuc()\n",
    "\n",
    "mf = scf.RHF(molecule)\n",
    "print(\"Running Restricted Hartree-Fock...\")\n",
    "mf.kernel()\n",
    "print(\"Running Restricted CCSD...\")\n",
    "mf_cc = cc.RCCSD(mf)\n",
    "mf_cc.kernel()\n",
    "print(\"Running Full Configuration Interaction...\")\n",
    "mf_fci = fci.FCI(mf)\n",
    "energy_FCI = mf_fci.kernel()[0]\n",
    "\n",
    "energy_RHF = mf.e_tot\n",
    "energy_RCCSD = mf_cc.e_tot\n",
    "\n",
    "print(\"Hartree-Fock energy: %20.16f Ha\" % (energy_RHF))\n",
    "print(\"CCSD energy: %20.16f Ha\" % (energy_RCCSD))\n",
    "print(\"FCI energy: %20.16f Ha\" % (energy_FCI))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be464515",
   "metadata": {},
   "source": [
    "所需要的单、双电子积分可以通过PySCF的ao2mo模块进行计算，详细内容可参考PySCF的Documentation相关章节："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "0620d102",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_orb = molecule.nao_nr()\n",
    "hcore = mf.get_hcore()\n",
    "mo_coeff = mf.mo_coeff\n",
    "one_body_mo = reduce(np.dot, (mo_coeff.T, hcore, mo_coeff))\n",
    "two_body_mo = ao2mo.restore(1, pyscf.ao2mo.get_mo_eri(\n",
    "    molecule, mo_coeff, compact=False), \n",
    "    n_orb\n",
    ")\n",
    "one_body_int = np.zeros([n_orb * 2] * 2)\n",
    "two_body_int = np.zeros([n_orb * 2] * 4)\n",
    "\n",
    "for p in range(n_orb):\n",
    "    for q in range(n_orb):\n",
    "        one_body_int[2 * p][2 * q] = one_body_mo[p][q]\n",
    "        one_body_int[2 * p + 1][2 * q + 1] = one_body_mo[p][q]\n",
    "        for r in range(n_orb):\n",
    "            for s in range(n_orb):\n",
    "                two_body_int[2 * p][2 * q][2 * r][2 * s] = two_body_mo[p][s][q][r]\n",
    "                two_body_int[2 * p + 1][2 * q + 1][2 * r + 1][2 * s + 1] = two_body_mo[p][s][q][r]\n",
    "                two_body_int[2 * p + 1][2 * q][2 * r][2 * s + 1] = two_body_mo[p][s][q][r]\n",
    "                two_body_int[2 * p][2 * q + 1][2 * r + 1][2 * s] = two_body_mo[p][s][q][r]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "得到单、双电子积分后，根据公式：\n",
    "\n",
    "$$\n",
    "\\hat{H} = \\sum_{p, q}{h^{p} _ {q} E^{p} _ {q}} + \\sum_{p, q, r, s}{\\frac{1}{2} g^{pq} _ {rs} E^{pq} _ {rs} }\n",
    "$$\n",
    "\n",
    "得到分子哈密顿量$\\hat{H}$。使用mindquantum的FermionOperator算符进行构造："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "309f3dd9",
   "metadata": {},
   "outputs": [],
   "source": [
    "hamiltonian_fermOp_1 = FermionOperator()\n",
    "hamiltonian_fermOp_2 = FermionOperator()\n",
    "\n",
    "for p in range(n_orb * 2):\n",
    "    for q in range(n_orb * 2):\n",
    "        hamiltonian_fermOp_1 += FermionOperator(\n",
    "            ((p, 1), (q, 0)),\n",
    "            one_body_int[p][q]\n",
    "        )\n",
    "for p in range(n_orb * 2):\n",
    "    for q in range(n_orb * 2):\n",
    "        for r in range(n_orb * 2):\n",
    "            for s in range(n_orb * 2):\n",
    "                hamiltonian_fermOp_2 += FermionOperator(\n",
    "                    ((p, 1), (q, 1), (r, 0), (s, 0)),\n",
    "                    two_body_int[p][q][r][s] * 0.5\n",
    "                )\n",
    "\n",
    "hamiltonian_fermOp_1 = normal_ordered(hamiltonian_fermOp_1)\n",
    "hamiltonian_fermOp_2 = normal_ordered(hamiltonian_fermOp_2)\n",
    "hamiltonian_fermOp = hamiltonian_fermOp_1 + hamiltonian_fermOp_2\n",
    "hamiltonian_fermOp += FermionOperator((), energy_nuc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "此处使用Jordan-Wigner变换将费米子算符映射为Qubit算符："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "02e48d5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "hamiltonian_QubitOp = Transform(hamiltonian_fermOp).jordan_wigner()\n",
    "n_qubits = count_qubits(hamiltonian_QubitOp)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69d3c231",
   "metadata": {},
   "source": [
    "耦合簇算符的构造如下所示。耦合簇算符，尤其是高阶激发算符，有多种不同的公式，对应的变分参数的个数、量子比特门的数量也不相同。此处采用下文代码中参考文献里的Spin-adapted激发算符，该公式也是[ADAPT-VQE](https://doi.org/10.1038/s41467-019-10988-2)的[开源代码](https://github.com/mayhallgroup/adapt-vqe)中所使用的："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "0b074fa7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def Pij(i: int, j: int):\n",
    "    ia = i * 2 + 0\n",
    "    ib = i * 2 + 1\n",
    "    ja = j * 2 + 0\n",
    "    jb = j * 2 + 1\n",
    "    term1 = FermionOperator(\n",
    "        ((ja, 0), (ib, 0)),\n",
    "        1.0\n",
    "    )\n",
    "    term2 = FermionOperator(\n",
    "        ((ia, 0), (jb, 0)),\n",
    "        1.0\n",
    "    )\n",
    "    return np.sqrt(0.5) * (term1 + term2)\n",
    "\n",
    "def Pij_dagger(i: int, j: int):\n",
    "    return hermitian_conjugated(Pij(i, j))\n",
    "\n",
    "def Qij_plus(i: int, j: int):\n",
    "    ia = i * 2 + 0\n",
    "    ib = i * 2 + 1\n",
    "    ja = j * 2 + 0\n",
    "    jb = j * 2 + 1\n",
    "    term = FermionOperator(\n",
    "        ((ja, 0), (ia, 0)),\n",
    "        1.0\n",
    "    )\n",
    "    return term\n",
    "\n",
    "def Qij_minus(i: int, j: int):\n",
    "    ia = i * 2 + 0\n",
    "    ib = i * 2 + 1\n",
    "    ja = j * 2 + 0\n",
    "    jb = j * 2 + 1\n",
    "    term = FermionOperator(\n",
    "        ((jb, 0), (ib, 0)),\n",
    "        1.0\n",
    "    )\n",
    "    return term\n",
    "\n",
    "def Qij_0(i: int, j: int):\n",
    "    ia = i * 2 + 0\n",
    "    ib = i * 2 + 1\n",
    "    ja = j * 2 + 0\n",
    "    jb = j * 2 + 1\n",
    "    term1 = FermionOperator(\n",
    "        ((ja, 0), (ib, 0)), \n",
    "        1.0\n",
    "    )\n",
    "    term2 = FermionOperator(\n",
    "        ((ia, 0), (jb, 0)), \n",
    "        1.0\n",
    "    )\n",
    "    return np.sqrt(0.5) * (term1 - term2)\n",
    "\n",
    "def Qij_vec(i: int, j: int):\n",
    "    return [Qij_plus(i, j), Qij_minus(i, j), Qij_0(i, j)]\n",
    "\n",
    "def Qij_vec_dagger(i: int, j: int):\n",
    "    return [hermitian_conjugated(i) for i in Qij_vec(i, j)]\n",
    "\n",
    "def Qij_vec_inner(a: int, b: int, i: int, j: int):\n",
    "    vec_dagger = Qij_vec_dagger(a, b)\n",
    "    vec = Qij_vec(i, j)\n",
    "    result = FermionOperator()\n",
    "    for i in range(len(vec)):\n",
    "        result += vec[i] * vec_dagger[i]\n",
    "    return result\n",
    "\n",
    "def spin_adapted_T1(i: int, j: int):\n",
    "    \"\"\"\n",
    "    Spin-adapted single excitation operators.\n",
    "    Reference:\n",
    "        Scuseria, G. E. et al., J. Chem. Phys. 89, 7382 (1988)\n",
    "    \"\"\"\n",
    "    ia = i * 2 + 0\n",
    "    ib = i * 2 + 1\n",
    "    ja = j * 2 + 0\n",
    "    jb = j * 2 + 1\n",
    "    term1 = FermionOperator(((ia, 1), (ja, 0)), 1.0)\n",
    "    term2 = FermionOperator(((ib, 1), (jb, 0)), 1.0)\n",
    "    return [term1 + term2]\n",
    "\n",
    "def spin_adapted_T2(creation_list: list, annihilation_list: list):\n",
    "    \"\"\"\n",
    "    Spin-adapted double excitation operators.\n",
    "    Reference:\n",
    "        Igor O. Sokolov et al., J. Chem. Phys. 152, 124107 (2020) \n",
    "        Ireneusz W. Bulik et al., J. Chem. Theory Comput. 11, 3171−3179 (2015)\n",
    "        Scuseria, G. E. et al., J. Chem. Phys. 89, 7382 (1988)\n",
    "    \"\"\"\n",
    "    p = creation_list[0]\n",
    "    r = annihilation_list[0]\n",
    "    q = creation_list[1]\n",
    "    s = annihilation_list[1]\n",
    "    tpqrs1 = Pij_dagger(p, q) * Pij(r, s)\n",
    "    tpqrs2 = Qij_vec_inner(p, q, r, s)\n",
    "    tpqrs_list = [tpqrs1, tpqrs2]\n",
    "    return tpqrs_list"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d336c7da",
   "metadata": {},
   "source": [
    "然后根据公式\n",
    "\n",
    "$$\n",
    "\\hat{T} = \\sum_{p\\not\\in occ\\\\\\\\q\\in occ}{\\theta^{p} _ {q} E^{p} _ {q}} + \\sum_{pq\\not\\in occ\\\\\\\\rs\\in occ}{\\theta^{pq} _ {rs} E^{pq} _ {rs}} + \\dots\n",
    "$$\n",
    "\n",
    "以及公式\n",
    "\n",
    "$$\n",
    "| \\Psi_{UCC} \\rangle = \\exp{(\\hat{T} - \\hat{T}^{\\dagger})} | \\Psi_{HF} \\rangle\n",
    "$$\n",
    "\n",
    "构造幺正耦合簇UCCSD所需要的单、双激发算符。同样地，最后使用Jordan-Wigner变换将构造好的费米子激发算符转化为Qubit算符："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "e0b95586",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_orb_occ = sum(molecule.nelec) // 2\n",
    "n_orb_vir = n_orb - n_orb_occ\n",
    "occ_indices = [i for i in range(n_orb_occ)]\n",
    "vir_indices = [i + n_orb_occ for i in range(n_orb_vir)]\n",
    "\n",
    "T1_singles = []\n",
    "T2_doubles = []\n",
    "constant = 1.0\n",
    "for p_idx in range(len(vir_indices)):\n",
    "    p = vir_indices[p_idx]\n",
    "    for q_idx in range(len(occ_indices)):\n",
    "        q = occ_indices[q_idx]\n",
    "\n",
    "        tpq_list = spin_adapted_T1(p, q)\n",
    "        for idx in range(len(tpq_list)):\n",
    "            tpq = tpq_list[idx]\n",
    "            tpq = tpq - hermitian_conjugated(tpq)\n",
    "            tpq = normal_ordered(tpq)\n",
    "            if len(list(tpq.terms)) > 0:\n",
    "                T1_singles.append(tpq)\n",
    "\n",
    "for p_idx in range(len(vir_indices)):\n",
    "    p = vir_indices[p_idx]\n",
    "    for q_idx in range(p_idx, len(vir_indices)):\n",
    "        q = vir_indices[q_idx]\n",
    "        for r_idx in range(len(occ_indices)):\n",
    "            r = occ_indices[r_idx]\n",
    "            for s_idx in range(r_idx, len(occ_indices)):\n",
    "                s = occ_indices[s_idx]\n",
    "\n",
    "                tpqrs_list = spin_adapted_T2([p, q], [r, s])\n",
    "                for idx in range(len(tpqrs_list)):\n",
    "                    tpqrs = tpqrs_list[idx]\n",
    "                    tpqrs = tpqrs - hermitian_conjugated(tpqrs)\n",
    "                    tpqrs = normal_ordered(tpqrs)\n",
    "                    if len(list(tpqrs.terms)) > 0:\n",
    "                        T2_doubles.append(tpqrs)\n",
    "\n",
    "uccsd_operator_pool_fermOp = T1_singles + T2_doubles\n",
    "uccsd_operator_pool_QubitOp = [Transform(op).jordan_wigner() for op in uccsd_operator_pool_fermOp]\n",
    "n_params = len(uccsd_operator_pool_fermOp)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8083669",
   "metadata": {},
   "source": [
    "Hartree-Fock初态波函数与前文类似，在Jordan-Wigner变换下可以通过$N$个$\\text{X}$门制备："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "4f423b71",
   "metadata": {},
   "outputs": [],
   "source": [
    "occ_indices_spin = [i for i in range(n_orb_occ * 2)]\n",
    "hartreefock_wfn_circuit = Circuit([X.on(i) for i in occ_indices_spin])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5836e4b2",
   "metadata": {},
   "source": [
    "幺正耦合簇波函数拟设对应的量子线路无法直接根据公式得到，需要作Trotter-Suzuki分解。Trotter-Suzuki分解可用下式描述：\n",
    "\n",
    "$$\n",
    "e^{\\hat{A} + \\hat{B}}\\approx(e^{\\hat{A}/k}e^{\\hat{B}/k})^{k}\n",
    "$$\n",
    "\n",
    "幺正耦合簇波函数拟设在该分解下可写成一系列指数算符的乘积：\n",
    "\n",
    "$$\n",
    "|\\Psi_{UCC} \\rangle=\\lim_{N\\rightarrow\\infty}{\\prod _ {k=1}^{N}{\\prod _ {i}^{N_{i}}{e^{\\frac{\\theta_{i}}{N}\\hat{\\tau} _ {i}}|\\Psi_{0} \\rangle}}}\n",
    "$$\n",
    "\n",
    "其中$\\hat{\\tau} _ {i}$属于反厄米算符$\\{ E^{p} _ {q} - E^{q} _ {p}, E^{pq} _ {rs} - E^{sr} _ {qp} \\}$之一。\n",
    "\n",
    "在Jordan-Wigner变换后，$\\hat{\\tau} _ {i}$被映射为一系列Pauli算符的直积的加和：\n",
    "\n",
    "$$\n",
    "\\hat{\\tau} _ {i} = \\sum_{\\mu}{c_{\\mu} P_{\\mu i}}\n",
    "$$\n",
    "$$\n",
    "P_{\\mu i} = \\prod^{\\otimes} _ {k}{\\sigma_{k \\mu i}}, \\sigma_{k \\mu i} \\in \\text{{X, Y, Z}}\n",
    "$$\n",
    "\n",
    "一串Pauli算符的直积的指数$\\exp{(P_{\\mu i})}$有等效的量子线路表示，此处调用mindquantum的`decompose_single_term_time_evolution`函数生成对应的线路："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "07f04eb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "ansatz_circuit = Circuit()\n",
    "ansatz_parameter_names = []\n",
    "for i in range(n_params):\n",
    "    qubitOp = uccsd_operator_pool_QubitOp[i]\n",
    "    ansatz_parameter_names.append(\"p\" + str(i))\n",
    "    for term, coeff in qubitOp.terms.items():\n",
    "        trottered_qubitOp = QubitOperator(term, coeff)\n",
    "        ansatz_circuit += decompose_single_term_time_evolution(\n",
    "            trottered_qubitOp,\n",
    "            {\"p\" + str(i): coeff.imag}\n",
    "        )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8647768b",
   "metadata": {},
   "source": [
    "和前文一样，总的波函数线路等于HF初态+UCC拟设："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "f5432480",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===============================Circuit Summary===============================\n",
      "|Total number of gates  : 16772.                                            |\n",
      "|Parameter gates        : 832.                                              |\n",
      "|with 44 parameters are : p41, p28, p38, p10, p13, p31, p16, p3, p4, p29... |\n",
      "|Number qubit of circuit: 12                                                |\n",
      "=============================================================================\n",
      "Number of parameters: 44\n"
     ]
    }
   ],
   "source": [
    "total_circuit = hartreefock_wfn_circuit + ansatz_circuit\n",
    "total_circuit.summary()\n",
    "print(\"Number of parameters: %d\" % (len(ansatz_parameter_names)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "需要注意的是，虽然此时变分参数的总数也是44个，但是单个变分参数对应的费米子算符和前文有所不同，故线路中总的量子比特门的数目也不一样。读者也可尝试不同的文献中给出的不同的激发算符来构造幺正耦合簇拟设所需的线路，以探索各种形式的激发算符对线路深度、变分参数数量的影响。使用新的`total_circuit`重新构造参数化量子线路算子`molecule_pqc`："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ae77444",
   "metadata": {},
   "outputs": [],
   "source": [
    "molecule_pqc = generate_pqc_operator(\n",
    "    [], ansatz_parameter_names, total_circuit,\n",
    "    Hamiltonian(hamiltonian_QubitOp))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86e28a29",
   "metadata": {},
   "source": [
    "此处我们定义一个函数`energy_objective`，该函数返回能量以及能量关于每个变分参数的导数。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "d86d8d6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def energy_objective(parameters: np.ndarray, \n",
    "                     pqc: mq.nn.pqc.PQC):\n",
    "    params_zero = ms.Tensor([[0]], ms.float32)\n",
    "    params_opt = ms.Tensor(parameters, ms.float32)\n",
    "    energy, _, grad = pqc(params_zero, params_opt)\n",
    "    return energy.asnumpy()[0, 0], grad.asnumpy()[0, 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56c7e48b",
   "metadata": {},
   "source": [
    "使用`energy_objective`得到初始能量。同样地，此时使用全为0的变分参数作为初始值："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "67ff6dd0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial energy:  -7.8633575439453125\n"
     ]
    }
   ],
   "source": [
    "init_amplitudes = np.zeros(n_params)\n",
    "initial_energy, grad = energy_objective(init_amplitudes, molecule_pqc)\n",
    "print(\"Initial energy: %20.16f\" % (initial_energy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最后，调用SciPy的minimize，使用BFGS优化器进行VQE优化模拟。不同的波函数拟设、不同的优化算法会对最终结果有一定的影响，但不会有过大的差别，一般会在小数点后三位以下。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "b1bf98fe",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Warning: Desired error not necessarily achieved due to precision loss.\n",
      "         Current function value: -7.882353\n",
      "         Iterations: 15\n",
      "         Function evaluations: 98\n",
      "         Gradient evaluations: 86\n",
      "Optimized energy:  -7.8823528289794922\n",
      "Optimized amplitudes: \n",
      " [-2.22631022e-04 -3.55662712e-02 -6.57172984e-16  3.97093689e-15\n",
      " -5.33872718e-15  4.49781015e-15  6.39930343e-04 -3.04572518e-03\n",
      "  1.90964239e-03  2.86849924e-04  1.60018334e-02  2.79161962e-14\n",
      " -8.25128341e-16  2.27834837e-16 -5.78023991e-16  9.87635287e-15\n",
      " -4.86620764e-16 -8.24096183e-16 -7.91120404e-16  8.74431219e-06\n",
      "  6.76903853e-04 -9.24311014e-04 -5.46154055e-02  9.22750289e-04\n",
      "  3.11477834e-03  1.42495172e-02 -2.14072217e-16  6.86482069e-18\n",
      "  4.65684505e-19 -2.58642139e-19 -1.66390756e-15  6.09039261e-17\n",
      "  2.24957290e-17  6.77074460e-16  9.03726787e-04  3.11791185e-03\n",
      "  1.42543419e-02 -1.55147586e-14  4.46844092e-16  1.97836927e-15\n",
      "  1.13411621e-15  4.25921631e-04  3.50493542e-04  5.40203928e-02]\n"
     ]
    }
   ],
   "source": [
    "result = minimize(energy_objective, \n",
    "                  init_amplitudes, \n",
    "                  args=(molecule_pqc),\n",
    "                  method=\"BFGS\",\n",
    "                  jac=True,\n",
    "                  options={\n",
    "                      \"disp\":True\n",
    "                  })\n",
    "print(\"Optimized energy: %20.16f\" % (result.fun))\n",
    "print(\"Optimized amplitudes: \\n\", result.x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7cd5d983",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
